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We present a systematic study of the influence of different forcing types on the statisti- 
cal properties of supersonic, isothermal turbulence in both the Lagrangian and Eulerian 
frameworks. We analyse a series of high-resolution, hydrodynamical grid simulations with 
Lagrangian tracer particles and examine the effects of solenoidal (divergence-free) and 
compressive (curl- free) forcing on structure functions, their scaling exponents, and the 
probability density functions of the gas density and velocity increments. Compressively 
driven simulations show a significantly larger density contrast, a more intermittent be- 
haviour, and larger fractal dimension of the most dissipative structures at the same root 
mean square Mach number. We show that the absolute values of Lagrangian and Eule- 
rian structure functions of all orders in the integral range are only a function of the root 
mean square Mach number, but independent of the forcing. With the assumption of a 
Gaussian distribution for the probability density function of the velocity increments on 
large scales, we derive a model that describes this behaviour. 
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1. Introduction 

Knowledge of the statistical characteristics of turbulence is a key prerequisite for un- 
derstanding turbulent flows on virtually all scales (Frisch 1995 Lesieur 19931. While 
common terrestrial flows are incompressible, astrophysical flows are highly supersonic 
and compressible. For example, the birth of stars in the interstellar medium is thought 
to be controlled by supersonic turbulence ( Mac Low fc Klessen|2004 Scalo & Elmegreen 
2004[|McKee fc Ostriker||2007[ ) . As turbulence is by definition a process characterised by 
a chaotic and irregularly fluctuating velocity field, there is a scale-dependent spatial and 
temporal correlation of fluid quantities (Ishihara et al. 2009). The scale-dependent en 



ergy, density and velocity distributions, and the dynamical evolution thereof, are strongly 
influenced by the properties of the turbulence. While large improvements were made in 
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the understanding of incompressible turbulence in the last few years (e.g. with new tech- 
niques like superstatistics and Lagrangian statistics, |Beck||2004| |Toschi fc Bodenschatz| 



2009 1, there are still open questions in our understanding of compressible turbulence. The 



non-local, inter-scale processes of compressible turbulence arising for example in shock 
fronts change the Richardson-Kolmogorov picture of the energy cascade of incompress- 
ible turbulence, where scale-locality is crucial for the existence of universal statistics in 
the inertial range. Basic questions as the existence of the inertial range, the associated 
scaling laws, and the influence of intermittency are still open. Guided by the idea of ba- 
sic physical quantities like the momentum and the kinetic energy, there is a trend using 
mass- weighted velocity increments to describe turbulence in a compressible medium. For 



example, |Kritsuk et al.\ ( |2007[ ), [Schmidt et al.\ ( |2008[ ), and [Galtier fc Banerjee] ( |2011 1 



describe the intermittency appearing in their supersonic, compressible, numerical simu- 
lations with a phenomenological model using mass-weighted velocity increments. Since 
Lagrangian tracer particles are designed to follow the turbulent flow and thus the mass 
flux in numerical simulations, they are a powerful tool to study compressible turbulence, 
with an intrinsic mass-weighting. Furthermore, phenomenological models treating inter- 
mittency arising in compressible turbulence were developed in the last years. Boldyrev] 



et al. (2002) used the multifractal phenomenological model of She & Leveque (1994), 



which describes the most dissipative structures in incompressible turbulence as filaments, 
and modified it such that the most dissipative structures are sheets instead of filaments. 
Nevertheless, despite of these improvements over the last years, the theoretical under- 
standing of compressible turbulence is still poor and analytical derivations are rare. 
Three spatial or temporal ranges have to be distinguished for turbulence: the viscous 
dissipation range at small scales, the inertial range at intermediate scales, and the range 
above the turbulent injection scale, the so-called integral range. The velocity structure 
functions, an average of the difference of two velocities separated by an increment in 
space or time, increase from the dissipation range to the integral range. Here, we also 
consider the statistics of turbulence in the integral range, where temporal correlations 
of the velocity field are exponentially damped, and the structure functions saturate over 
several integral time scales. A complete understanding of the statistical properties in this 
range is necessary as a solid foundation for further studies of the complex and compli- 
cated behaviour of turbulence in the inertial range. In this paper we shed some light 
on this unattended topic. We compute Lagrangian statistics of the density and velocity 
fields of supersonic turbulence, using data from high-resolution, three-dimensional simu- 
lations of driven turbulence with up to 1024'^ grid cells and up to 512'^ tracer particles. 
By calculating probability density functions and structure functions, we analyse the sen- 
sitivity of compressible turbulence on the forcing of the turbulence. In particular, we 
compare two limiting cases of purely solenoidal (divergence-free) and purely compressive 
(curl-free) forcing. We present a detailed analysis of the scaling properties of the velocity 
structure functions and intermittency in the inertial subrange. We show that a simple 
analytic formula describes the saturated structure functions of all orders in the integral 
range, with the r.m.s. Mach number as single input parameter. This formula can be used 
to check the convergence of the structure functions. 

In §[2] we explain the numerical setup, describe the implementation of the different forc- 
ings, the tracer particles, and define the structure functions and the statistical moments 
used to analyse the simulations. We analyse structure functions, their scaling proper- 
ties, intermittency, and the probability density function (PDF) of the mass density and 
the velocity increments in § |3] In § |4j we present a simple analytic formula describing 
the behaviour of the saturated structure functions of all orders in the integral range. A 
summary of our results and conclusions are given in § |5] 
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2. Simulation and Methods 

We solve the hydrodynamical equations on a uniform grid with 256^, 512'^, and 1024'^ 
grid points, using the piecewise parabolic method (PPM, Colella & Woodward 1984), 
implemented in the grid code FLASH3 ( |Fryxell et al]|2000[ (Pubey et aL|2008[ ). We start 



with gas of uniform density at rest and uniformly distributed tracer particles also at rest. 
We place one tracer particle in every other grid cell, such that the simulations contain 
128"^, 256'^, and 512^ tracer particles, respectively. Since we assume isothermal gas, the 
pressure, P = pcs^, is proportional to the density p with the fixed sound speed Cg. We 
solved the continuity equation and the Euler equation with a stochastic forcing term F, 

^%(vV)s--V-v, (2.1) 



dt 



at 



+ (v • V)v = -c^Vs - 



(2.2) 



where v is the velocity field and s = ln(/9/(p)v') is the natural logarithm of the mass 
density divided by the mean (volume-weighted) mass density. The parameters of the 
simulation are: {p)v = 1, Cg = 1 and the computational domain has a box length L = 1 
with periodic boundary conditions. The numerical simulation is evolved for ten dynam- 
ical time scales T — L/2V where V is the integral velocity and the relevant quantities 
are stored in intervals of 0.01 T for the tracer particles and in intervals of 0.1 T for the 
grid. We do not apply physical viscosity, so we have to rely on purely numerical viscos- 
ity. This requires resolution studies. Hence we investigate our results for different grid 
resolutions: 256^, 512^, and 1024^ (Appendix |a| . It can be shown that the numerical 
viscosity of PPM can be used as an implicit way of treating physical viscosity, as long as 
a large-enough scale separation is guaranteed (Benzi et aL||2008|. 



2.1. Forcing module 

The random forcing term F is derived from a stochastic Ornstein-Uhlenbeck process 



with finite autocorrelation time scale, T^c (Eswaran & Pope 1988 Schmidt et al. 2009 



[Federrat h et aL|2010 ). It gives a stochastic force field F that varies smoothly in space and 
time. The Ornstein-Uhlenbeck process generates the forcing in Fourier space (fc-space) 
by solving a differential equation. 



dF{k,t) = Fo{k,T,c)rH^) 



dW(t) 



F(k,t) 



d^ 
T 



(2.3) 



where the dW(t) is a three-dimensional Gaussian random increment with zero mean and 
standard deviation dt, generated by a Wiener process. 'P''(k) is a projection tensor in 
Fourier space. In index notation, this operator is 

7'«„-(k) - CV^.,{k) + (1 - OVhjik) , (2.4) 

where "P^ ~ Sij — kikj/k^ and 7^" — kikj/k^ are the fully solenoidal and compressive 
projection operators respectively. By setting ( — 1, the forcing field is purely solenoidal 
(i.e., V • F = 0). In contrast, setting C = 0, the forcing field is purely compressive (i.e., 
V X F = 0). The natural mixture of forcing modes is obtained for ( — 0.5 (Federrath 
2010), which leads to a velocity distribution of (vfi) / (v'^ot) = 1/3- To investigate 



et al. 



the iiinuence of the different forcings, we focus on the limiting cases of purely solenoidal 
forcing (C — 1) and purely compressive forcing (C = 0). 

The forcing amplitude -Fb(k, T) is a three-dimensional parabolic function, only containing 
the large (integral) scales 1 <| k |< 3, peaking at fc = 2, which corresponds to half of the 
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box size L/2, as we measure k in units of 2t:/L. The amplitude of the forcing is adjusted, 
so that in both cases the volume-weighted r.m.s Mach number is My ~ 5.5, when the 
state of stationary, fully-developed turbulence is reached. 



The last term in equation (2.3) is a stochastic diffusion term that ensures exponential 
decrease of the autocorrelation function of the forcing. We set the autocorrelation time 
Tac of the forcing equal to the dynamical time scale T. 

The density and velocity statistics of turbulence produced by a mainly compressive force 



field are investigated in Schmidt et al. ( 2009, ), while a systematic statistical comparison of 



solenoidal and compressive forcing is discussed in Schmidt et al. { 2008 1 and in Federrath 



et a/.| ( |2008||2009[|2010 ). 

2.2. Tracer particles 

We start with uniformly distributed tracer particles at rest. Afterwards they can move 
freely within the computational domain. The velocity and density of the tracer particles 
are calculated with a cloud-in-cell interpolation from the grid at the beginning of each 
time step. Given the interpolated velocity, the tracer particles are then moved with 
an Euler method, based on the hydrodynamical time step. The tracers thus follow the 
gas flow in the Lagrangian frame of reference. Instead of the linear interpolation of 
the neighbouring grid points, a second-order (triangular-shaped cloud) and third-order 
(tricubic) space interpolation, as well as a higher-order temporal integration scheme 
(predictor-corrector type) were tested, but they did not lead to statistically significant 
differences. As the tracer particles have no influence on the fluid, they are passive tracers 
of the fluid motion. 

2.3. Velocity increments and structure functions 

In order to calculate the increments, we use the following definition of the time-dependent, 
Lagrangian velocity increment 

svr{t,T) = vrit+T)-vr{t), (2.5) 

where r is a temporal increment and vY^{t) is the velocity in spatial direction i € {x, y, z} 
of the mth tracer particle at the time t. The space-dependent, Eulerian velocity incre- 
ments are defined as 

Svr{r,e) = vr{r + e)~v^{r) (2.6) 

Sv^"'{r,£)^vj;^{r + £)-v]^{r), (2.7) 

where r is the position of the tracer n, £ is the spatial increment between the tracer 
particles m and n, and wy = v ■ £ with £ — £/£ being the unit vector in the direction £. 
The Lagrangian structure function (LSF) 

LSP{T)^{{\6v:\t,T)nn^ + {\Svl-{t,T)nn.+ {\SvT{t,T)n^)t/i (2.8) 

is obtained by averaging the velocity increments over the different tracer particles to, the 
three directions of the coordinate system x, y, z and over t e [2, 10] T. This is reasonable 
because of the time invariance and isotropy in the state of fully-developed turbulence. 
In practise, we randomly select 5 x 10^ tracer particles for all 801 time samples with 



t ^ 2T in the fully-developed state for the averaging procedure (2.8 1. We checked the 
validity of this approach by doing these calculations also with all 512'^ tracer particles 
for one time-line starting at t = 2r to ensure that the used number of sampling pairs 
has no statistically significant influence on our results (see Appendix |b]) . For calculating 
the Eulerian structure functions (ESFs) 

ES^{1) ^ ((|(5z;™"(r,^)r)^„ + ([^^^^^"(r, £)|f + (|<5<"(r,^)|P)^„),/3 (2.9) 
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ESP{£)= {\Svl^^{r,£)\P)„,,,,t (2.10) 

the simulation box was divided in 16"^ sub boxes. For m, a fixed number of tracer particles 
is chosen homogeneously distributed over all sub boxes. To obtain a constant sampling of 
the ESF with £, for each m, a subset oc 1 /r^ of tracer particles of every sub box is selected 
for n, where r is the distance from m to the centre of the sub box. As the number of sub 
boxes increases proportional to r^, this procedure ensures that for each lag £, roughly 
the same number of sampling pairs is used. The selection procedure is normalised in a 
way that for m, nearly the same number of tracer particles is selected as for n. The ESFs 
are calculated for 81 snapshots in time intervals of At — T/10 each with about 10^° 
sampling pairs. We tested our results with different numbers of sub boxes, where with 
insufficient sub boxes (< 8"^) we have to calculate the ESFs with many more sampling 
pairs to get a good statistic on the smallest scales. Using more than 16'^ sub boxes showed 
no effective improvement of the distribution. With IG'^ sub boxes, we calculate the ESFs 
with different numbers of sampling pairs to ensure that increasing the sampling pairs has 
no statistically significant influence on our results. We provide detailed convergence tests 
in Appendix [B] Since all increments are calculated on the tracer particles, the structure 
functions are intrinsically mass-weighted. 

2.4. Statistical Moments 

In order to calculate the higher-order moments of the PDFs in § [3| we use the following 
definition for the first four standardised central moments: 

mean: {q)=^qp{q)Aq (2.11a) 

standard deviation: dq = \J {{q — (q))^) (2.116) 

skewness: S„ = 1^^''^'' (2.11c) 

kurtosis: /C„ = ^^^^^^^ , (2.11d) 
a 

where Aq is the bin width of the PDF p{q). With this definition a Gaussian has a skewness 
Sq — and a kurtosis ICq = 3. 



3. Results 



As discussed in |Federrath et a/.| p009| |2010[ ) and [Price fc Federrath| ( |2010[ ), the fluid 

reaches an equilibrium state of fully-developed, supersonic turbulence after about two 
turbulent crossing times, t k, 2T. We therefore restrict our analysis to times t ^ 2T. 
Figure [T] shows the time evolution of the mass-weighted r.m.s Mach number and the 
mass-weighted averaged density, calculated with the data of the tracer particles. As the 
tracer particles are advected by the flow, their density is correlated with the mass density. 
As a consequence, the average of a quantity over all tracer particles is mass-weighted. We 
denote quantities calculated with this quasi Lagrangian statistics with a subscript M. 
Figure[l]demonstrates that at t « 2 T, a regime of statistically fully-developed turbulence 
is reached. The right panel of figure [T] shows that the compressive forcing yields a nearly 
eight times larger mass- weighted mean density and fluctuations thereof. The mean values 
indicated in figure [T] are averaged over the interval 2T ^t ^ 10 T, and the errors are the 
1 a variations in time. In the state of fully- developed turbulence, the time-averaged r.m.s 
Mach number on the tracer particles is (A^soi,*/)* = 4.9 (subscript t for time average) 
with solenoidal forcing, and (A^comp,M)t = 4.4 with compressive forcing. In Federrath 
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Figure 1. Mass-weighted r.m.s Mach number (left), and averaged mass- weighted density (right), 
as a function of the dynamical time, calculated by averaging over all tracer particles. Within 
the first 2T, a. statistically steady state was reached for both solenoidal (sol) and compressive 
(comp) forcing. The mean values were averaged over t £ [2, 10] T and the errors are the 1 a 
variations in time. 



\et aL|p008[ table 1), the volume- weighted Mach number of this simulation was measured 
{{■Msoi,v)t = 5.3, {Mcomp,v)t = 5-6). Clearly, the intrinsic mass-weighting of the tracer 
particles influences the statistical properties, discussed in more detail below. The values of 
the volume- weighted and mass- weighted r.m.s. Mach number for the different forcings are 
summarised in table [T] (first row). To illustrate the different flow patterns for solenoidal 
and compressive forcing, we show a randomly selected slice through the mid plane of 
the computational domain at t = 6 T in figure [2j It shows the logarithm of the mass 
density computed on the grid cells in this slice, as well as the norm of the velocity of the 
tracer particles that are in one slice with thickness of 0.1 grid cells. Each dot represents 
a tracer particle with the colour corresponding to the norm of the velocity. The density 
fluctuations are more space-filling with solenoidal forcing, and have smaller amplitude, 
while compressive forcing yields larger voids and denser regions. The bottom panel of 
figure [2] shows a magnification of a head-on collision of two flows that leads to a strong 
shock in the simulation with compressive forcing. The density field shows a sharp, well- 
defined shock front. In these compressed regions, the tracer particles accumulate and 
have a significantly lower Mach number, Mm ~ 1 — 2. These stagnation points, showing 
a correlation of high density with low velocities ( Ballesteros-Paredes et a/. [2003 Klessen 
et al. 2005 e.g.), are important in the theory of star formation ( Mac Low & Klessen 2004 ), 



as they are good candidates for the formation of dense-cores, which are the progenitors 
of individual stars and binary stellar systems. This correlation in the stagnation points 
causes the mass-weighted values of the Mach number to be smaller than the volume- 
weighted ones. Compressive forcing excites more head-on collisions and shock fronts, so 
this effect has a stronger influence in that case. 



3.1. The probability density function of the gas density 
The probability density function (PDF) of the gas density p{p) and its standard devia- 



tion (T^ are important quantities in astrophysics. For instance, |Padoan & Nordlund|(|2002|) 



and Hennebelle & Chabrier ( 2008[ 2009) relate the density PDF to the mass distribu 
tion of dense gas cores and stars. Padoan et al. (1997) and Passot & Vazquez-Semadeni 



( 1998 ) have shown that the standard deviation grows proportional to the Mach number 



of the turbulent flow, if the density PDF is close to a log-normal distribution (see Price 



et aL||20lT for a recent, extended study). Federrath et al. (20101 demonstrated that the 
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Figure 2. Top panels: Slice in the xy-plane at z = 0.5 and time t = 6T as an example for 
the regime of statistically fully-developed, compressible turbulence. The logarithm of the mass 
density of the grid cells in this slice (blue) as well as the norm of the velocity of the tracer 
particles (red) in one slice with thickness of 0.1 grid cell are displayed for solenoidal (left) and 
compressive forcing (right). Bottom right: Magnified slice of a shock front in the simulation with 
compressive forcing (see white box in upper right panel). 



density PDF is not only influenced by the r.m.s Mach number, but also by the forcing 
parameter C, and presented a modification of the existing expression, which takes the 
ratio of solenoidal and compressive modes of the forcing into account. In many numerical 
experiments of driven, supersonic, isothermal turbulence with solenoidal and/or weakly 
compressive forcing, it was found that the density PDF is close to a log-normal distri- 
bution (e.g. , Padoan et a/.||1997 KlessenpOOO Lemaster fc Stone]|2008 Federrath et al. 



2008), 



p{s) = 



27r(T, 



exp 



2a? 



(3.1) 
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where s = hi[p/ {p)v) is the logarithm of the density divided by the volume-weighted 
mean density. Li et al.\ ( 2003 ) argued that the mass- weighted density distribution is also a 



log-normal, with the same standard deviation as the volume-weighted distribution. With 
the assumption of a log-normal density PDF, the authors derived a relation between the 
mass-weighted and the volume-weighted quantities. 



IM 



(3.2) 



Figure [3] presents the time-averaged {t e [2, 10] T) PDF of the quantity sm, calculated on 
the tracer particles in linear (left panel) and logarithmic scaling (right panel). The first 
four standardised central moments (see section 2.4) of the density PDF pm{s) are sum- 
marised in table [T] together with the volume- weighted moments calculated in [Federratli 
et al. (2010 table 1). Compressive forcing yields a significantly broader mass- weighted 
density distribution with standard deviation, (Ts,m about 1.5 times larger at the same 
volume-weighted Mach number than solenoidal forcing ((Ts,soi,Af — 1-23, comp,Af = 
1.77). The peak is shifted to larger values of the logarithmic density ((s)soi,j\/ — 0.81, 
(s)comp,M — 2.37). The PDF is compatible with a Gaussian distribution for solenoidal 
forcing. However, the Gaussian fit (black dashed line) shows that the density PDF has 
weak non-Gaussian contributions in the wings of the distribution. On the other hand, 
in the density PDF obtained from compressive forcing, the discrepancy to the Gaus- 
sian distribution in both wings is more prominent. The deviations from the log-normal 
distribution for the compressive forcing is caused by both a physical and a numerical 



effect. Price & Federrath (2010) analysed the influence of measuring the density PDF 



by using a grid-based simulation and an SPH simulation and found that the PDF of the 
SPH particles increases slightly in the high-density tail with increasing resolution and 
decreases in the low-density tail. We expect that this effect will decrease the deviations 
from a log-normal distribution for both forcing types in our simulation as the resolution is 



increased. Federrath et al. (2010 figure 6) analysed the volume-weighted density PDFs 
for resolutions of 256'*, 512'^, and 1024^, showing that changing the resolution affects 



solenoidal and compressive forcing in roughly the same way. Schmidt et al. (2009) also 



argue that the deviations from a log-normal PDF produced by compressive forcing are 
a genuine effect. From this we can conclude that the stronger non log-normal features 
seen for compressive forcing likely have a physical origin rather than a purely numerical 
one. The higher moments of the distribution with compressive forcing (iScomp,Af = —0.57, 
^comp,A/ — 3.50) show larger deviations from the Gaussian values (5 = 0, /C = 3) than for 
solenoidal forcing {SsoI.m — —0.13, /Csoi.Af = 2.95). Checking the relation (3.2) between 



the mean value and the standard deviation demonstrates that for solenoidal forcing, 
the assumption of a log-normal PDF is nearly fulfilled (cr^ ^/2 — 0.76). In contrast, 
we find larger discrepancy for compressive forcing {(^1 ^omv mI'^ — 1-57). Measuring the 



volume- weighted PDF, Federrath et al. (2010) also reported small deviations from a 
Gaussian distribution for solenoidal forcing (5soi,y = —0.10, /Csoi,v = 3.03) and larger 
deviations for compressively-driven turbulence (5comp,y = —0.26, /Ccomp.y = 2.91). The 
mass-weighted quantities show larger discrepancy from the Gaussian values than the 
volume- weighted quantities. 



Padoan et al. (1997) and Passot & Vazquez-Semadeni (1998) motivated a linear relation 



between the r.m.s Mach number and the standard deviation of the linear density. 



<^p/Po = bM, 



(3.3) 



with a proportionality constant b. Several measurements of b were obtained from different 
simulations and with different types of forcing, and yielded b between 0.26 and 1.05 
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M 



solenoidal compressive 
mass-weighted volume-weighted mass-weighted volume-weighted 



4.9 ±0.2 
0.81 ±0.04 
1.23 ±0.03 
-0.13 ±0.05 
2.95 ±0.07 



5.3 ±0.2 
-0.83 ±0.05 

1.32 ±0.06 
-0.10 ±0.11 

3.03 ±0.17 



4.4 ±0.3 
2.37 ±0.24 
1.77 ±0.09 
-0.57 ±0.17 
3.50 ±0.33 



5.6 ±0.3 
-3.40 ±0.43 

3.04 ± 0.24 
-0.26 ± 0.20 

2.91 ± 0.43 



Table 1. Mass- weighted and volume-weighted r.m.s. Mach number (first row) and the first 
four standardised central moments of the PDF Pm(s) (column 1 and 3) and Pv{s) (column 2 
and 4) for solenoidal and compressive forcing. These quantities are averaged in time when the 
equilibrium state of fully-developed, supersonic turbulence is reached (see i|3]and j ]3.1[ ) and the 
errors are the standard deviation in time. 




Figure 3. Mass-weighted PDF Pm{s) of the logarithmic mass density s = ln(p/(p)v) for 
solenoidal (solid line) and compressive (dashed-dotted line) forcing in linear (left panel) and 
logarithmic scaling (right panel), calculated on the tracer particles. The PDFs are calculated 
for 81 time steps in the state of fully-developed turbulence t ^ 2T and averaged. The error bars 
indicate the standard deviation of the temporal fluctuations. The dashed lines show log-normal 
flts with the mean value and the standard deviation as fitting parameters. 



(Passot & Vazquez-Seinadeni||1998 Li et 112003 Kritsuk et a/. 112007 Beetz et a/. 112008 



Federrath et al. 2008 Schmidt et al. 2009 Price et al. 2011). From the distribution 



Pm{s) shown in figure |3j we calculate pm{p)- Since equation (3.3) was derived with 



volume-weighted quantities, we have to transform our results using the relation (see, Li 
eFaI][2003| , 

, , dMdy , , ,^ ^, 

PMiP) cc (X ppv[p) , (3.4) 

where pm (p) and pv (p) are the mass- weighted and the volume- weighted PDFs of the gas 
density, respectively. We find Up^soiy = 1-90 for purely solenoidal and Cpxomp.v — 6.03 for 
purely compressive forcing. With the volume-weighted Mach number of this simulation, 
we get 6soi = 0.36 for solenoidal forcing and &comp = 1-08 for compressive forcing, in good 



agreement with Federrath et al. (2008). 



3.2. The probability density functions of velocity increments 

The simplest set of correlation functions to quantify the statistical properties in a com- 
pressible, supersonic turbulent flow is the distribution of the velocity increments and its 



higher moments, the structure functions, defined by (2.8 to 2.101. The deviation of the 
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Figure 4. Top panel: PDFs of the velocity increment Svi in the Lagrangian framework for 
solenoidal (left panel) and compressive forcing (right panel). It shows the PDFs with five different 
temporal increments r G {0.01, 0.08, 0.4, 2, 4}T. Bottom panel: PDFs of the velocity increment 
Svi in the Eulerian framework for solenoidal (left panel) and compressive forcing (right panel). It 
shows the PDFs with six different spatial increments £ £ {0.006, 0.02, 0.06, 0.12, 0.25, 0.49}L. 
For small temporal or spatial lags, t and £, respectively, the PDFs differ from a Gaussian distri- 
bution because of intermittency. For large r or £, they converge towards a Gaussian distribution. 



structure function scaling exponents from the predicted values of the Kolmogorov model 
( Kolmogorov||194l" ) is an effect of intermittency (e.g., She fc Leveque|[l994 ) . A property 
of intermittency is that the PDFs of the velocity fluctuations become more and more 
non-Gaussian on smaller and smaller scales ( Gotoh et aL||200^ Mordant et aL|[2"002 ). 
Figure [4] shows the PDF of the velocity increment Svi in the Lagrangian framework for 
five temporal increments r e (0.01, 0.08, 0.4, 2, 4}T and in the Eulerian framework for 
six spatial increments £ € {0.006, 0.02, 0.06, 0.12, 0.25, 0.49}L. The PDFs follow a Gaus- 
sian distribution for r — ^ T and £ L. Decreasing the spatial or temporal increment, 
the Gaussian PDFs vary continuously towards distributions with exponential wings in- 
dicating the intermittent behaviour of the turbulent velocity field. Figure [5] shows the 
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Figure 5. Kurtosis as a function of the temporal or spatial increment r (left) and I (right). 
The solid and dashed-do tted lines are respectively, the values for solenoidal and compressive 
forcing, calculated from (3.51. The crosses and stars are respectively the values for solenoidal 
and compressive forcing, calculated with the five PDFs shown in figure|4] The horizontal, dotted 
line is the value of a Gaussian distribution. 



kurtosis (see § 2.4) of the distributions of the Lagrangian (left panel) and Eulerian (right 
panel) velocity increments, calculated with the structure functions 



K{T)^LS\T)l{LS\T)f 



(3.5) 



(solid and dash-dotted lines) and the values calculated with the PDFs (crosses and stars) . 
The kurtosis can be used as a measure for the deviations of the distributions of the 
velocity increments from a Gaussian distribution. In the Lagrangian framework, the 
kurtosis obtained with solenoidal forcing converges towards the Gaussian value (/C — 3) 
on times comparable with the dynamical time scale, t k, \ T . The compressive forcing 
yields PDFs converging already on smaller temporal lags, r « 0.7 T than solenoidal 
forcing. Compressive forcing develops a more intermittent behaviour with larger kurtosis 
than the solenoidal forcing for times r < 0.08 T. As the non-Gaussian wings of PDFs 



of the density and/or velocity is caused by intermittency (see Federrath et al. 2010), 



this analysis of the kurtosis and the more intermittent behaviour of the compressive 
forcing confirm our observation of the density PDF and its deviation from the log-normal 
distribution (see 



3.1) 



In the Eulerian framework, compressive forcing yields a more intermittent behaviour with 
a larger kurtosis on nearly all spatial scales than solenoidal forcing. However, the kurtosis 
obtained with compressive forcing converges at the same scale {£ « 0.23 L) towards the 
Gaussian value as the kurtosis of the solenoidal forcing. Comparing the kurtosis of the 
Lagrangian and Eulerian structure functions strengthens the conclusion that Lagrangian 
statistics are more intermittent than Eulerian ones (here shown for two limiting types of 



forcing) as already observed by Benzi et al. (2010), but only for purely solenoidal forcing 



In the region, where the kurtosis is converged towards the Gaussian value (r ^ 2. ST for 
the Lagrangian framework and £ > 0.4 L for the Eulerian framework), we average the 
structure functions to calculate the mean value of the saturated structure functions, as 
discussed in § 13.31 
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3.3. Lagrangian and Eulerian structure functions 

Figure [6] shows the LSF and ESF up to order p = 7 for solenoidal and compressive 
forcing. We calculate the saturation values of the structure functions on the largest 
scales by averaging them in the range t G [2.5, 5] T for the LSFs and £ £ [0.4, 0.7] L 
for the ESFs. The result is displayed as black lines in figure [6) The compressive forcing 
yields structure functions that converge to Gaussian values already on smaller scales, as 
observed in figure [5] For the two forcings, the saturation values are different. This can be 
explained with the different mass- weighted r.m.s Mach number A4 m, observed in figure[l| 
because the structure functions S of order p for infinite increments are 5*^(00) oc A^^^ 
(see discussion in §|4]). 

In turbulence theory, the structure functions follow a power law in the inertial range 



(3.6) 



with the scaling exponents ^{p) and C(p)- To calculate these scaling exponents, we use 



the inertial range as constrained by [Federrath et a^.| pOOOt |2010| (0.067 < i/L < 0.2), 
and transform it to the Lagrangian framework (0.16 < t/T < 0.34) with r cx i"^/^. This 
relation follows directly from the Kolmogorov four-fifth law (e.g., [FHsch||1995 ), implying 
that the third-order structure function scales linearly with i in the Eulerian framework. 
In Burgers turbulence, t (x £, which follows from the assumption of a constant averaged 
energy transport through the scales, e oc E{£)/t (x v{tf' /t, and assuming Burgers scaling 
for the second-order velocity increment, 5v'^{t) oc v'^{t) oc I. Using Burgers scaling for t in 
the transformation of the inertial range leads to (0.067 < t/T < 0.2) in the Lagrangian 
framework. The so-called method of extended self-similarity (ESS) proposed by 'Benzi 
et al. ( 1993 1 allows for an increased scaling range between the smallest scales, influenced 



by the resolution, and the largest scales with a direct influence of the forcing. Using ESS, 
we thus extend the fitting range to 0.067 < t/T < 0.34, which covers both the trans- 
formation with T oc i'^^^ and t (x £. For the Eulerian structure functions, we extended 
the scaling range to (0.05 £/L < 0.22) for which we obtain a reasonable power-law 
scaling with ESS. Figure [T] shows the ESS scaling plots, i.e., plots of the logarithm of 
the structure functions calculated with equations (2.8) and (2.9) for the different orders 
as a function of the logarithm of the second- and third-order structure function in the 
Lagrangian and Eulerian framework, respectively. The black lines indicate linear fits for 
the ESS-measurement of the relative scaling exponents, 



Zl{p) 



m. 



Ze{p) = 



CM 

C(3) 



(3.7) 



which are summarised in table [2] for the Lagrangian framework (2th and 3th columns) 
and the Eulerian framework (5th and 6th columns) for solenoidal and compressive forc- 
ing, respectively. To compare our results with data from incompressible turbulence, we 
refer in tabl e [2| to the data of numerical simulations of subsonic turbulence published by 
Benzi et al.\ ([2010| table 2, ReynoldsA ~ 600, BBFLTIO) and|Gotoh eFaI](p002] table 3, 



ReynoldsA = 381, GFN02). We only show the data of the transverse structure functions 
of GFN02 in the Eulerian framework, because the differences between the longitudinal 
and transverse structure function are negligible compared to the differences between the 
supersonic and the subsonic results. We expect that our results of the structure function, 
averaged over the three directions of the coordinate system are in-between the results of 
the longitudinal and transverse structure function. Additionally, we compare our results 
calculated with tracer particles, discussed here, with the results of the same simulations, 
but calculated with the p^^^ mass-weighted velocities measured on the grid published by 



Lagrangian and Eulerian analysis 



13 







Lagrangian 






Eulerian 






p 


ZLip) 


Zl{p) 


BBFLTIO 


Ze(p) 


Zsip) 


SFK08 


SFK08 


GFN02 




sol. 


comp. 




sol. 


comp. 


sol. 


comp. 


(transversal) 


1 


0.57 ±0.02 


0.55 ±0.04 




0.422 ± 0.008 


0.45 ±0.02 


0.539 


0.605 


0.369 ± 0.004 


2 


1.0 


1.0 


1.0 


0.75 ±0.02 


0.78 ±0.03 


0.840 


0.869 


0.701 ±0.01 


3 


1.30 ±0.05 


1.30 ±0.09 




1.0 


1.0 


1.0 


1.0 


0.998 ±0.02 


4 


1.51 ±0.07 


1.5 ±0.1 


1.66 ±0.02 


1.20 ±0.04 


1.13 ±0.07 


1.080 


1.066 


1.26 ±0.03 


5 


1.67 ±0.09 


1.6 ±0.1 




1.36 ±0.05 


1.19 ±0.09 


1.112 


1.100 


1.49 ±0.04 


6 


1.8 ±0.1 


1.6 ±0.2 


2.10 ±0.10 


1.49 ±0.07 


1.2 ±0.1 






1.69 ±0.05 


7 


1.9 ±0.1 


1.6 ±0.2 




1.61 ±0.09 


1.2 ±0.2 






1.86 ±0.05 



Table 2. Relative sc aling exponents o f the structure functions LS''{t) and ES^{t) calculated 
with the ESS method (Benzi et al. 1993 1 for solenoidal and compressive forcing in the Lagrangian 
(2th and 3th columns) and Eulerian (5th and 6th columns) frameworks. For a comparison with 
the scaling exponents of the structure functions in an incompressible turbulent medium, we refer 



to t he data of numerical simulations published by Benzi et a/!| (2010l (BBFLTIO, 4th column) 
and |Gotoh et a l.\ (|2002| (GFN02, 9th column). The 7th and 8th columns show the relative 
scaling exponents of the transversal structure function of our simulation, but calculated with 
the ve locities measured on the grid and mass-weighted with p^^^ published in 
(2008t(SFK08). 



Schmidt et al. 



Schmidt et al. (2008) (SFK08, 7th and 8th columns). These scaling exponents are also 
calculated for the transverse structure functions. The mass- weighting p^/^, used by many 



authors (see e.g.,|Kritsuk et a/.|2007 


Kowal & Lazarian 2007, 


Schmidt et al. 2008 


Galtier| 


& Banerjee|2011 ), follows from the assumption of a constant 


mean volume energy trans- 



we can analyse the influence of the different forcings in each framework. The relative 
scaling exponents show only a significant difference between the scaling behaviour of the 
solenoidal and compressible forcing in the Lagrangian framework for the highest order. 
The scaling exponents of the compressive forcing are slightly below the scaling exponents 
of the solenoidal forcing for higher orders. However, in the Eulerian framework, this effect 
is stronger and the compressive forcing causes scaling exponents to stay nearly constant 
above an order p > i, such that there is a significant difference between the scaling expo- 
nents of the solenoidal and compressive forcing. With the measured scaling exponents, we 
can quantify the intermittency in the different frameworks by calculating the differences 



to the predicted Kolmogorov (1941 K41) scaling. In the Lagrangian framework and for 
the highest order, p = 7, the scaling exponents are 46 ± 3% and 54 ± 6% smaller than 
the K41 value for the solenoidal and compressive forcing, respectively. In the Eulerian 
framework, the scaling exponents are 31 ± 4% and 49 ± 9% smaller than the K41 value. 
For solenoidal forcing, the scaling exponents show a more intermittent behaviour in the 
Lagrangian framework than in the Eulerian one. This is consistent with our analysis of 
the kurtosis in figure [s] and the results of Benzi et al. (20101. For compressive forcing, 



we have an intermittency of the same order for both frameworks. The intense density 
fluctuations in the simulation with compressive forcing cause a more intermittent be- 
haviour and scaling exponents that deviate stronger from the K41 values. The stronger 
influence of the compressive forcing on the intermittency in the Eulerian framework is 
an important result that needs further studies. 

To estimate the influence of shocks and other non-local, inter-scale processes arising in su- 
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0.1 1.0 0.1 

T / T I / L 

Figure 6. Lagrangian (left) and Eulerian (right) structure functions up to the order of p = 7 for 
solenoidal and compressive forcing, calculated with the velocity increment of the tracer particles. 
The error bars of the structure functions indicate the standard deviation in time. To calculate 
the absolute values of the saturated structure functions, the Lagrangian structure functions were 
averaged in the range r £ [2.5, 5] T, and the Eulerian structure functions were averaged in the 
range I £ [0.4,0.7] L, indicated by the horizontal black lines in the integral range. 



personic, compressible turbulence, we compare our results with the data of other subsonic, 
incompressible simulations. Our scaling exponents are below those for incompressible me- 
dia in both frameworks and are significantly different from the data of BBFLTIO and 
GFN02. This indicates a more intermittent behaviour in our supersonic, compressible, 
turbulent flow, which is even stronger for the simulations with compressive forcing. The 
comparison of the scaling exponents of the tracer particles with the p^^^ mass-weighted 
results of the grid in the Eulerian framework shows that the p^/^ multiplier does not have 
the same effect as the averaging over tracer particles. The intrinsic mass-weighting of the 
tracer particles shows a less intermittent behaviour for both forcings than the results of 
SFK08. 



3.4. Intermittecy models for inertial range scaling 

With the relative scaling exponents of the last section, we can compare our results with 
the predictions of intermittency models. We use the generalised equation of the phe- 



nomenological model of She & Leveque (1994) introduced by Dubrulle (1994) for the 



scaling exponents in the Eulerian framework 
Zs(p)-(1-Ab)| + 



A, 



(1 



(3.8) 



With the assumptions r oc (."^/^ and LS^ oc (eT^^)TP/^, where (e^) are the moments of the 
energy dissipation at the time scale r, one can show a similar equation for the Lagrangian 



framework, using the same arguments and derivations of She & Leveque (1994): 



(l-/3f) 



(3.9) 



For simplicity we use r cx p/'^ for the transformation into the Lagrangian framework, 
instead of r oc £, and treat the influence of compressibility in both frameworks by having 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 

logio(LS'(T)) log,o(ES'(e)) 

Figure 7. Extended self-similarity for the Lagrangian (left panel) and Eulerian (right panel) 
structure function for solenoidal and compressive forcing. The black lines indicate a linear fit 
for the calculation of the relative scaling exponents, which are summarised in table [2] 



different values for A and /3 compared to the K41 theory (see e.g., Boldyrev et aL|[2002 



Schmidt et al. 2008). Figure [s] shows the measured scahng exponents and the fits with 



equation (3.8) and (3.9). For the fitting procedure we follow the idea of Schmidt et al. 
(2008 2009) and set A = 1, which follows from Burgers scaling, r cx as used in 
the last section, leaving us with only one free fitting parameter. With the measured /? 
we can calculate the co-dimension of the most dissipative structures C = A/(l — /3), 
which is connected to the actual dimension of the most dissipative structures via D = 
3 — C. The latter quantifies how volume-filling the most dissipative structures are in the 
turbulent medium. From our fits we get Dl.soI — 0.87, -Di.comp = 1-17 in the Lagrangian 
framework and sol = l-H, £'£;,comp = 1-55 in the Eulerian one. In the Eulerian 
framework, the most dissipative structures are between filamentary structures {D = 1, 
as in|She & Leveque (1994)) and sheet-Hke structures {D = 2, as proposed by Boldyrev 



et al. \ (2002 1 for the Kolmogorov-Burgers model). In the Lagrangian framework and for 
solenoidal forcing, the most dissipative structures are close to filamentary structures. 
Compressive forcing yields a larger fractal dimension than solenoidal forcing in both 
frameworks. Although the whole turbulent flow is more space filling for solenoidal forcing, 
as observed in figure [2] the most dissipative structures of the compressive forcing have a 
larger dimension and are thus more space filling. However, it is unclear how to interpret 
these results in the one-dimensional Lagrangian framework of temporal increments rather 
than spatial increments as in the Eulerian framework. In the Eulerian framework, we can 
compare our results with the dimensions we get by calculating the scaling exponents with 
the mass- weighted velocities of the grid. Schmidt et al. (2008) measured Dgoi = 1-82 and 
-Dcomp = 1.92, showing the same trend between the solenoidal and compressive forcing, 
but larger than the values we measured on the tracer particles. The reason for these 
differences is the more intermittent behaviour of the scaling exponents, as discussed 
above (see also table [2| . 



16 



L. Konstandin, C. Federrath, R. Klessen and W. Schmidt 



■3 

M 



2.5 
2.0 

1 .5 

1 .0 

0.5 







Lagrangian 




)K SOL 

fit SOL 

COMP 
- fit COMP 



Eulerian 




4 

P 



4 
P 



Figure 8. Relative scaling exponents for the Lagrangian (left panel) and Eulerian (right panel) 
structure function for solenoida l and compressiv e forcing. The lines ind icate a fit with the in - 
termittency model proposed by Dubrulle (19941 with the assumption of Schmidt et al. (20081, 
A = 1. In the Lagrangian framework, we get a dimension of the most dissipative structures 
Dl,eoI ~ 0.87 and DL.comp = 1-17 for solenoidal and compressive forcing respectively. In the 
Eulerian framework we get De,soI ~ 1-11 and DE,comp = 1.55. The compressive forcing and 
the associated stronger density fluctuations cause a higher dimension of the most dissipative 
structures. 



4. A statistical theory of the large-scale velocity increments 

In this section, we show that the statistical properties of the velocity increments in a 
turbulent flow on large scales can be described with only one parameter, the r.m.s Mach 
number. This is valid for velocity increments in the Lagrangian and Eulerian framework. 



The structure functions defined by (2.8 1 to (2.101 can be expressed as the moments of 
the PDFs of the velocity increment, which are functions of t or £, so we can write for a 
general structure function, 

SP{a)^ I \6vf P{5v.a)d{5v), (4.1) 



where P{Sv, a) is the probability density of Sv with the increment a. In the last sec- 
tion, we showed that the PDFs of the velocity increments converge towards a Gaussian 
distribution on the largest scales. The Gaussian form can be understood analytically as 
a consequence of the central limit theorem, assuming that the two velocities, w'"(r + i) 
and u"(r) in space or v'^{t + r) and v"^{t) in time, are independent for large spatial or 
temporal increments. 

With the Gaussian assumption, we can express the structure functions on large scales as 

oo 

SP{a ^ cx)) = [ {Svfe-'i^ d {6v) (4.2a) 



r(2±i 





1 

2 



{V2a)P, (4.26) 



where SP{a) stands for any structure function of (2.81 to (2.101, a is the temporal or 
spatial increment, a is the standard deviation of the Gaussian distribution, and F is the 
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Gamma function. Equation ( |4.2p ) describes the moments of the Rayleigh distribution, 
which is also the result for the moments of the total structure function with a velocity 

increment Sv = \l &v'^ + 5vy + 5v1, if the increments 5vi follow a Gaussian distribution. 



Stutzki et al.\ ( |1998[ ) showed that 

{{5v{i~^^)f)^2Ml,cl, (4.3) 

where they used homogeneity and the fact that the autocorrelation vanishes for large 
spatial increments. In our case, the quantity M.m is a mass-weighted value, because the 



average in (4.3) is taken over the velocity increments of the tracer particles. Furthermore, 
we assume that the second-order structure function is proportional to the kinetic energy 
for large increments and as the longitudinal structure function and the structure function 
averaged over the three directions of the coordinate system have only one-third degree 
of freedom compared with the total structure function, 

{{6vf) = {6vl + 5vl + 6vl) = i{{5v,f) = 3((J«||)2) . (4.4) 



If we combine this with (4.2 ) and (4.3 1, we get a relation between the standard deviations 
of the Gaussian distributions and the r.m.s Mach number Mm'- 

{{5vf) = 2Ml,cl = 3^2 = 3^2 . (4.5) 



The second-order moment can thus be used as a normalisation for our formula ( |4.2[ ) to 
predict the saturation level of the pth-order structure function 



, 2 



S^{a^^)^^^{-^MMr . (4.6) 

VTT V3 



Figure [9] shows the structure functions of figure[6] but renormalised with equation (4.6 1 to 



the r.m.s. Mach number of the solenoidal forcing. The differences between the structure 
functions, driven by solenoidal and compressive forcing, vanishes in the integral range, 
what implies that the different forcings have no influence on the statistical properties 
of the structure functions in the integral range. Additionally, we verify this model by 
calculating the saturation behaviour with the measured r.m.s Mach number and compare 
the result with the saturation values extracted from figure [6j The result is summarised 
in figure |10| The measurements show an excellent agreement with the predicted values, 
for both solenoidal and compressive forcing. 



5. Summary and conclusions 

We have investigated the influence of solenoidal (divergence-free) and compressive 
(curl-free) forcing on the structure functions and density PDFs of a supersonic, com- 
pressible, turbulent flow using tracer particles in a set of three-dimensional numerical 
simulations. We analysed the density PDF, the PDFs of velocity increments, and the 
structure functions in the Lagrangian and Eulerian frameworks. As all of these quanti- 
ties were measured on tracer particles, we analysed mass-weighted statistics. Our main 
results and conclusions are: 

• The solenoidal forcing yields a density PDF close to a log-normal distribution. In con- 
trast, the compressive forcing yields distributions of the mass density that show stronger 
deviations from the log-normal shape in the wings of the distribution. 

• The compressive forcing excites stronger head-on collisions and shock fronts, which 
show a correlation between high density and low velocity, affecting the mass-weighted 
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Figure 9. Same as figure |6j but with structure functions of tlie compressive forcing normalised 
to tiie r.m.s. Macfi number of the solenoidal forcing, using our model prediction (4.61. The 
differences of the structure functions between the solenoidal and compressive forcing thus vanish 
in the integral range. 
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Figure 10. Values of the saturated structure functions me asure d from the simulation (crosses, 
stars, and error bars) and the values predicted from formula (4.6 1 (solid and dashed-dotted lines) 
for different orders p. The values are shown to the left of their respective order and the crosses 
and stars are from left to right the values of the Lagrangian, Eulerian, and the longitudinal 
Eulerian structure function. The error bars of the measured saturation values are the averaged 
errors of the structure functions. 



r.m.s Mach number, sucti that it becomes smaller than the volume-weighted Mach num- 
ber. Similar holds for the solenoidal forcing, but the effect is weaker, as solenoidal forcing 
yields smaller density contrasts at the same r.m.s. Mach number. 

• The Lagrangian framework exhibits a more intermittent behaviour than the Eulerian 
one, measured with the deviations of the relative scaling exponents from the predicted 
intermittency-free K41 values and also with the kurtosis as an example for the higher 
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moments of the PDF of the velocity increments. This analysis also shows that the tur- 
bulent medium, driven by the compressive forcing, is more intermittent than a medium, 
driven by solenoidal forcing. A comparison with simulations of incompressible turbulence 
shows that intermittency is stronger in a supersonic, compressible medium. 

• The influence of the different forcings are stronger in the Eulerian than in the La- 
grangian framework, measured with the relative scaling exponents and the kurtosis of 
the velocity increments. 

• The fractal dimension of the most dissipativc structures are larger for compressive 
than for solenoidal forcing. The dimensions are around those of filamentary structures 
(D = 1) in the Lagrangian framework and between filamentary and sheet-like structures 
(1 < _D < 2) in the Eulerian framework. Although the whole turbulent flow driven by the 
solenoidal forcing is more space filling the most dissipativc structures of the compressive 
forcing are more space filling. 

• The behaviour of the structure functions of all orders on the largest scales is deter- 
mined by the r.m.s. Mach number of the system. With the assumption of a Gaussian 
distribution for the velocity increments on the largest scale, we derived an analytic re- 
lation, predicting the absolute values of the structure functions in the integral range 
of turbulence. The statistical properties of the velocity increments are indistinguishable 



in the integral range for both forcings types and follow our derived formula (4.6 1. The 
different values of the saturated structure functions observed in figure [lO] are caused by 
the different Mach numbers {A4so\,m ~ 4.9 ± 0.2 and 7Mcomp,A/ — 4.4 ± 0.3), and not 
by different statistical properties obtained by the different forcings. Thus, the predic- 



tions based on equation (4.6) are independent of the energy injection mechanism of the 



turbulence, but only depend on the mass-weighted r.m.s. Mach-number. 
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Appendix A. Influence of the numerical grid resolution 

Figure [II] (left panel) shows the LSF of order p = 7 for both forcing types and with the 
grid resolutions 256^, 512^, and 1024^. The right panel shows the same for the ESF. The 
LSFs are calculated with 128'^ « 2.1, 512^^ « 16.8 and 5 million tracer particles for the 
different grid resolutions, respectively. The ESFs are calculated with 16^ sub boxes and 
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Figure 11. Structure functions of order p = 7 for different grid resolutions, shown for both 
forcing types. The LSF (left panel) with a low grid resolution of 256^ and 512'^ were calculated 
with 128^ and 256'^ tracer particles and the simulation with 1024^ grid cells was calculated with 
5 X 10^ tracer particles. The ESF (right panel) was calculated with 16"^ sub boxes and 10^" 
sampling pairs. The structure functions of the compressive forcing was multiplied with a factor 
of 10 so that the structure functions are distinguishable. 



with 10^" sampling pairs. Figure 11 shows that the structure functions of order p = 7 
differ by about 15%, caused by the different grid resolutions. This is of the same order 
as the 1 a variations in time of the structure functions indicated as error bars in figure [6] 
Therefore, the influence of the resolution is smaller than the temporal variations. 



Appendix B. Convergence test for the structure functions 

In order to verify that our calculations arc converged with a sufficient number of data 
pairs to sample the structure functions, we show that the structure functions do not 
change significantly by further increasing the number of sampling pairs. As large velocity 
fluctuations have a stronger influence on the higher orders of the structure functions 
and these events are very rare, the statistical convergence of the higher orders is slower 
compared with the lower orders. Thus, if we can demonstrate convergence for the higher 
order structure functions, this automatically holds for all lower orders. Figure 12 (left 
panel) shows the Eulerian structure function of order p — 7 for solenoidal and compressive 
forcing. The structure function of the compressive forcing is multiplied by a factor of 10, 
so that the structure functions of the different forcings are distinguishable. In order 
to check the convergence, we use one random time sample {t = AT) in the state of 
fully-developed turbulence, 16^ sub boxes and different numbers of sampling pairs (10^, 
10"'^°, 10"'^^). Increasing the number of sampling pairs further only influences small scales, 
£ < 0.07 i. The structure functions are converged on larger scales. For the Eulerian 
structure function, we also verified that the method of selecting tracer particles for the 
calculation with our procedure of sub boxes has no significant infiuence on the results. 
Therefore, we calculated the structure function with 10^^ sampling pairs and different 



numbers of sub boxes (8'^, 16^ and 32^). Figure 12 (right panel) shows that further 
increasing the number of sub boxes also only infiuences scales i < 0.05 L, and with 16'^ 
sub boxes, the structure functions are converged. 
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Figure 12. Left panel: ESF of order p — 7 for both forcing types, 16^ sub boxes and different 
numbers of sampling pairs. The structure functions are converged on large scales for more 
than 10^° samphng pairs. Right panel: Same as left panel, however with 10^^ sampHng pairs 
and different numbers of sub boxes. In both figures, the structure functions of the compressive 
forcing was multiplied by a factor of 10. 

For the Lagrangian structure function, we also have to verify that the structure functions 
do not change significantly with the number of sampling pairs. We calculate the LSF for 
all 512^ tracer particles for one time-line from t — 2T to t = 10 T and compare it 
with the LSF calculated with 5 and 10 million tracer particles. The results are shown in 
figure |13[ where the structure functions with compressive forcing are multiplied with a 
factor of 10. The structure functions calculated for solenoidal forcing are converged on 
all scales, and the structure functions calculated for compressive forcing show only small 
variations with the number of sampling pairs. The reason for the large fluctuations in the 
integral range in figure [13] is that the LSF was here calculated with one time-line only. 
Figure [13] shows that the time evolution of the forcing module has a direct infiuence on 
the amplitudes of the velocity increments in the integral range, but these fluctuations 
are smaller than the variations in time, we use as errors in figure J6] However, this direct 
influence vanishes on average by using different staring times for calculating the LSF. In 
the inertial range with t < 1 T, the structure functions in figure ]6] have a factor of about 
700 more sampling pairs for each bin. This large statistic we used there ensures that our 
structure functions are also converged in the inertial range. 
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factor of 10. 
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